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Abstract. The kinetics of nonequilibrium Bose-Einstein condensates are consid- 
ered within the framework of the Gross-Pitaevskii equation. A systematic derivation 
is given for weak small-scale perturbations of a steady confined condensate state. 
This approach combines a wavepacket WKB description with the weak turbulence 
theory. The WKB theory derived in this paper describes the effect of the conden- 
sate on the short-wave excitations which appears to be different from a simple 
renormalization of the confining potential suggested in previous literature. 



1 Introduction 

Bose-Einstein condensate (BEC) was first observed in 1995 in atomic vapors 
of 87 Rb [1], 7 Li [2] and 23 Na [3]. Typically, the gas of atoms is confined by 
a magnetic trap [1], and cooled by laser and evaporative means. Although 
the basic theory for the condensation was known from the classical works of 
Bose [4] and Einstein [5] , the experiments on BEC stimulated new theoretical 
work in the field (an excellent review of this material is given in [6]). 

A lot of theoretical results about condensate dynamics are based on the 
assumption that the condensate band can be characterized by some temper- 
ature T and chemical potential /i, the quantities which are clearly defined 
only for gases in thermodynamic equilibrium. Often, however, the condensa- 
tion is so rapid that the gas is in a very nonequilibrium state and hence, one 
requires the use of a kinetic rather than a thermodynamic theory [9-11]. An 
approach using the quantum kinetic equation was developed by Gardiner et 
al [9,10] who used some phenomenological assumptions about the scattering 
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amplitudes. Phenomenology is unavoidable in the general case due to an ex- 
treme dynamical complexity of quantum gases the atoms in which interact 
among themselves and exhibit wave-particle dualism. Most phcnomenologi- 
cal assumptions are intuitive or arise from a physical analogy and are hard 
to validate (or to prove wrong) theoretically. In particular, it was proposed 
that the ground BEC states act onto the higher levels via an effective po- 
tential. In the present paper we are going to examine this assumption in a 
special case of large occupation numbers, i.e. when the system is more like a 
collection of interacting waves rather than particles and which allows a sys- 
tematic theoretical treatment. In what follows we show systematically that 
such an assumption is not true for such systems. For dilute gases, with a 
large number of atoms at low temperatures, one obtains the Gross-Pitaevskii 
(GP) equation for the condensate order parameter [7,8]: 

idti) + AV - |V| V - Uip = 0, (1) 

where the potential U is a given function of coordinate, see for example figure 
1. We emphasize that the area of validity of GP equation is restricted to a 
narrow class of the low-temperature BEC growth experiments and the latest 
stages in other BEC experiments. However, we will study the GP equation 
because it provides an important limiting case for which one can rigorously 
test the phenomenological assumptions made for more general systems. We 
would like to abandon the approach where the system is artificially divided 
into a T = condensate state and a thermal "cloud" because this "cloud" in 
reality is far from the thermodynamic equilibrium and we believe that this 
fact affects the BEC dynaimcs in an essential way. As in many other non- 
equilibrium and turbulent systems, fluxes of the conserved quantities through 
the phase space are more relevant for the theory here than the temperature 
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and the chemical potential. Performance of a thermodynamic theory here 
would be as poor as a description of waterfalls by a theory developed for 
lakes. 1 Again, the GP equation is used in our work for both the ground and 
the excited states which limits our analysis only to the low temperature and 
high occupation number situations. 

In fact the idea of using GP equation for describing BEC kinetics is not 
new and it goes back to work of Kagan et al [11], who used a kinetic equation 
for waves systematically derived from the GP equation ignoring the trap- 
ping potential and assuming turbulence to be spatially homogeneous [12]. A 
similar method has been used to investigate optical turbulence [13]. Clas- 
sical weak turbulence theory yields a closed kinetic equation for the long 
time behavior of the energy spectrum without having to make unjustifiable 
assumptions about the statistics of the processes [14-16,18,22,24]. Second, 
the kinetic equation admits classes of exact equilibrium solutions [14,19,20]. 
These can be identified as pure Kolmogorov spectra [12-14], namely equilib- 
ria for which there is a constant spectral flux of one of the invariants, the 
energy, 



A very important property of the particle cascade is that it transfers the 
particles to the small k values (inverse cascade). This transfer will lead to an 
accumulation at small fc's which is precisely the mechanism of the BE con- 
densation, see figure 1. The energy cascade is toward high values of k which 

eventually will lead to "spilling" over the potential barrier corresponding to 

1 This comparison was suggested by Vladimir Zakharov to illustrate irrelevance of 
the thermodynamic approach to the turbulence of dispersive waves. 




and the "number of particles" , 
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an evaporative cooling, see figure 1. After the formation of strong condensate 
one can no longer use weak turbulence theory, as the weak turbulence theory 
assumes small amplitudes. However, one can reformulate the theory using a 
linearization around the condensate, (as oppose to linearization around the 
state), as in [13]. Consequently this changes the dominant system interactions 
from 4-wave to 3-wave processes. 

Kolmogorov-type energy distributions over the levels (scales) are dramat- 
ically different from any thermodynamic equilibrium distributions. Thus, the 
condensation and the cooling rates will also be significantly different from 
those obtained from theories based on the assumptions of a thermodynamic 
equilibrium and the existence of a Boltzmann distribution. As an example, 
a finite-time condensation was predicted by Kagan, Svistunov and Shlyap- 
nikov [11], whose work was based on the theory of weak homogeneous turbu- 
lence. 

However, application of the theory of homogeneous turbulence to the 
GP equation has its limitations. Indeed, when the external potential is not 
ignored in the GP equation, the turbulence is trapped and is, therefore, in- 
trinsically inhomogeneous (e.g. a turbulent spot). Additional inhomogeneity 
of the turbulence arises because of the condensate, which in the GP equation 
case is itself coordinate dependent. This means, in particular, that the the- 
ory of homogeneous turbulence cannot describe the ground state effect onto 




Fig. 1. Turbulent cascades of energy E and particle number N. 
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the confining properties of the gas and thereby test the effective potential 
approach. The present paper is aimed at removing this pitfall via deriving an 
inhomogeneous weak turbulence theory. 

The effects of the coordinate dependent potential and condensate can 
most easily be understood using a wavepacket (WKB) formalism that is ap- 
plicable if the wavepacket wavelength I is much shorter than the characteristic 
width of the potential well L, 

£ = y < 1. 

The coordinate dependent potential and the condensate distort the wavepack- 
ets so that their wavenumbers change. This has a dramatic effect on nonlinear 
resonant wave interactions because now waves can only be in resonance for 
a finite time. The goal of our paper is to use the ideas developed for the 
GP equation without the trapping potential and to combine them with the 
WKB formalism in order to derive a weak turbulence theory for a large set 
of random waves described by the GP equation. 

Note that idea to combine the kinetic equation with WKB to describe 
weakly nonlinear dynamics of wave (or quantum) excitations is quite old and 
can be traced back to Khalatnikov's theory of Bose gas (1952) and Landau's 
theory of the Fermi fluids (1956), see e.g. in [27]. It has also been widely used 
to describe kinetics of waves in plasmas, e.g. [28-31]. For plasmas, such a 
formalism was usually derived from the first principles. However, only phe- 
nomenological models based on an experimentally measured dispersion curves 
have been proposed so far for the supcrfluid kinetics. In this paper, we offer 
for the first time a consistent derivation starting from the GP equation which 
allows us to correct the existing BEC phenomenology at least for the special 
cases when the GP equation is applicable. 
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Technically, the most nontrivial new element of our theory appears through 
the linear dynamics (WKB) whereas modifications of the nonlinear part (the 
collision integral) are fairly straightforward. Thus, we start with a detailed 
consideration of the linear dynamics in section 2. Previously, linear excita- 
tions to the ground state were considered by Fetter [17] who used a test 
function approach to derive an approximate dispersion relation for these ex- 
citations. Fetter pointed out an uncertainty of the boundary conditions to be 
used at the ground state reflection surface. The WKB theory for BEC which 
is for the first time developed in the present paper allows an asymptotically 
rigorous approach which, among other things, allows to clarify the role of 
the ground state reflection surface. Indeed, as we will sec in section 3, the 
WKB theory is essentially different in the case when the condensate ground 
state is weak and can be neglected from the case of strongly nonlinear ground 
state. No suitable WKB description exists for the intermediate case in which 
the linear and the nonlinear effects are of the same order. However, in the 
Thomas-Fermi regime the layer of the intermediate condensate amplitudes 
is extremely narrow due to the exponential decay of the amplitude beyond 
the ground state reflection surface. This allowed us to combine the two WKB 
descriptions into one by formally re-writing the equations in such a way that 
they are correct in the limits of both weak and strong condensate. These 
equations will be wrong in the thin layer of intermediate condensate ampli- 
tudes, but this will not have any effect on the overall dynamics of wavepackcts 
because they pass this layer too quickly to be affected by it. 

In section 4 for the first time we present a Hamiltonian formulation of the 
WKB equations and derive a cannonical Hamiltonian the form of which is 
general for all WKB systems and not only BEC. The Hamiltonian formulation 
is needed to prepare the scene for the weak turbulence theory. In section 5 
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we apply weak turbulence theory to write a closed kinetic equation for wave 
action. This kinetic equation has a coordinate dependence of the frequency 
delta functions. Notice that coordinate dependence of the wave frequency has 
a profound effect on the nonlinear dynamics. The resonant wave interactions 
can now take place only over a limited range of wave trajectories which makes 
such interactions similar to the collision of discrete particles. 

2 Linear dynamics of the GP equation 

We will now develop a WKB theory for small-scale wave-packets, described 
by a linearized GP equation, with and without the presence of a background 
condensate. As is traditional with any WKB-type method we assume the 
existence of a scale separation e <C 1, as explained in section 1. In this 
analysis we will take I ~ 1 so that any spatial derivatives of a given large- 
scale quantity (e.g. the potential U or the condensate) are of order e. The 
transition to WKB phase-space is achieved through the application of the 
Gabor transform [23], 

g(x,k,t) = f f(e*\x-x \)e ik < x - x ^g(x ,t)dx , (2) 

where / is an arbitrary function fastly decaying at infinity. For our purposes 
it will be sufficient to consider a Gaussian of the form 

where d is the number of space dimensions. The parameter e* is small and 
such that Hence, our kernel / varies at the intermediate-scale. A 

Gabor transform can therefore be thought of as a localized Fourier transform, 
and in the limit e* — > becomes an exact Fourier transform. Physically, 
one can view a Gabor transform as a wavepacket distribution function over 
positions x and wavevectors k. 



8 Lvov, Nazarenko, West 

2.1 Linear theory without a condensate 

Linearizing the GP equation, to investigate the behavior of wavepackets tp 
without the presence of a condensate, we obtain the usual linear Schrodingcr 
equation: 

id t tp + AV> - Ui/> = 0, (3) 

where U is a slowly varying potential. Let us apply the Gabor transformation 
to (3). Note that the Gabor transformation commutes with the Laplacian, so 
that AW = AW. Also note that 

UW ~ UW + i(V x U)Vk&, 

where we have neglected the quadratic and higher order terms in e because 
W changes on a much shorter scale than the large scale function U. Combin- 
ing the Gabor transformed equation with its complex conjugate we find the 
following WKB transport equation, 

AM 2 = o, (4) 

where 

D t = d t + x-V + ie-d k , 

represents the total time derivative along the wavepacket trajectories in 
phase-space. The ray equations are used to describe wavepacket trajectories 
in (k,x) phase-space, 

x = dku, k = —Vuo. (5) 

The frequency uj, in this case, is given by u = k 2 + U, (again we use the 
notation k = |fc|). Equations (4) and (5) are nothing more than the famous 
Ehrenfest theorem from quantum mechanics. According to (5), the wavepack- 
ets will get reflected by the potential at points r# where U{rji) = k^^. We 
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will now move on to consider linear wavepackets in the presence of a back- 
ground condensate. 

2.2 Wavepacket dynamics on a condensate background 

One of the common assumptions in the BEC theory is that the presence of a 
condensate acts on the higher levels by just modifying the confining potential 
U, see for example [25]. If this was the case, the linear dynamics would still 
be described by the Ehrenfest theorem with some new effective potential. We 
will show below that this is not the case. 

Let us define the condensate ip as a nonlinear coordinate dependent so- 
lution of equation (1), with a lengthscale of the order of the ground state size 
(although it does not need to be exactly the same as the ground state). In 
what follows, we will use Madelung's amplitude-phase representation for ij) , 
namely 

Wa = VpW e i9 , (6) 

where v = 2V0 is the macroscopic speed of the condensate. It is well known 
that in this representation p obeys a continuity equation, 

Pt + div(pv) = 0. (7) 

For future reference, one should note that the second term in this expression 
is 0(e 2 ). Thus, p t is 0(e 2 ) too and it must be neglected in the WKB theory 
which takes into account only linear in e terms. We start by considering a 
small perturbation ^<1, such that 

^ = ^o(l + # (8) 

Substituting (8) into (1) we find 

id t <j> + A<j) + 2^ • - gU + </>* + 2|0| 2 + cf 2 + \<t>\ 2 <t>) = 0. (9) 
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where g = g(x) = \ipo\ 2 is a slowly varying condensate density. 

In a similar manner to the previous subsection, the rest of this derivation 
consists of Gabor transforming (9), combining the result with its complex 
conjugate and finding a suitable waveaction variable such that the transport 
equation represents a conservation equation along the rays. Such a deriva- 
tion is given in Appendix A. It yields to the following expression for the 
waveaction, 

(10) 

where 5ft and 3 mean the real and imaginary parts respectively. As usual, the 
transport equation takes the form of a conservation equation for waveaction 
along the rays, 

D t n(x, k, t) = 0, (11) 

where 

D t = d t + x- V + k-d k , (12) 
is the time derivative along trajectories 

x = dkUJ, k = -Vw. (13) 

The frequency is given by the following expression, 

lu = k^/k 2 + 2g. (14) 

One can immediately recognize in (14) the Bogolubov's formula [21] which 
was derived before for systems with a coordinate independent condensate 
and without a trapping potential. It is remarkable that presence of the po- 
tential U does not affect the frequency so that expression (14) remains the 
same. Obviously, the dynamics in this case cannot be reduced to the Ehren- 
fest theorem with any shape of potential U. Therefore, an approach that 



n(k, x, t) 



1 ojp 
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5ft0- 
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models a condensate's effect by introducing a renormalized potential would 
be misleading in this case. 

3 Applicability of WKB descriptions 

In this section we will investigate the applicability of the above theory. Let us 
consider a condensate which is a solution of the eigenvalue problem dttpo = 
—ifiipQ. Therefore, the GP equation (1) becomes 

QiPo + A?/>o - Q^o - UtPo = 0. (15) 

3.1 Weak condensate case 

Firstly, let us consider the case of a weak condensate so that the effect of 
the nonlinear term is small in comparison to the linear ones, \gipa\ < |A^o|- 
Since Q is a constant we observe that the Laplacian term acts to balance 
the external potential term (like in the linear Schrodinger equation) and the 
nonlinear term can be, at most, as big as the linear ones 

n „ 1 „ U(r ) > g, 

where r is the characteristic size of the condensate (it is defined as the 
condensate "reflection" point via the condition Q = U(ro), see below). 




Fig. 2. Regions of Applicability of WKB Descriptions. 



12 Lvov, Nazarenko, West 

Now for a WKB description to be valid we require kr n » 1, i.e. we require 
the characteristic length-scale of our wavepackets to be a lot smaller than that 
of the large-scales. Using this fact we find 

e » L „ u(r ) > g. 

Therefore, the condensate correction to the frequency, given by (14), is small. 
In other words the wavepacket docs not "feel" the condensate. Indeed, from 
hmax = U(rn) we have U{rn) 3> U(r ) and this implies that r« » r (where 
is the wavepacket reflection point, see figure 2). Thus, the condensate in 
this case occupies a tiny space at the bottom of the potential well and hence 
does not affect a wavepacket's motion. Therefore, a wavepacket moves as a 
"classical" particle described by the Ehrenfest equations (4) and (5). In fact, 
in this case it would be incorrect to try to describe the small condensate 
corrections via our WKB approach because these corrections are of order 
g ~ e 2 (the e 2 terms being ignored in a WKB description) . 

3.2 Strong condensate case 

Now we will consider a strong condensate such that 

f2^U + g»^ (16) 
I Vol 

i.e. the r dependence of the potential U is now balanced by the nonlinearity. 
This is usually referred to as the Thomas- Fermi limit [6]. Wavepackets now 
"feel" the presence of a strong condensate if g <~ k 2 . We see that the WKB 
approach is applicable because 

h 2 n ^ 1 |AVo| 



I Vol ' 

According to the ray equations a; is a constant along a wavepacket's trajec- 



tory, so we can find the packet's wavenumber from k 2 = \J g 2 + uo 2 — g. One 
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can see that k 2 remains positive for any value of g which means that the 
presence of the condensate does not lead to any new wavepacket reflection 
points (i.e. when k takes a value of zero). Thus, turbulence is allowed to 
penetrate into the center of the potential well. However, the group velocity 
increases when the condensate becomes stronger, d^ui ~ yj~p. This means that 
the density of wavepackets decreases toward the center of well. Therefore, the 
condensate tends to push the turbulence away from the center, toward the 
edges of the potential trap. 

To summarize, in the presence of a strong condensate we have two regions 
of applicability for our WKB descriptions, see figure 2. Wavepackets at a 
position r < r , in the central region of the potential well will evolve according 
to the WKB-condensate description (10) - (14). The Laplacian term only 
becomes important for r > tq where g is exponentially small. In this case the 
Ehrenfest description is appropriate. It will be shown in the next section that 
these two WKB descriptions can be combined into a single set of formulae. 

3.3 Unified WKB description 

It is interesting that taking the limit of zero condensate amplitude in the 
waveaction (10) results in the waveaction \\^\ 2 of the Ehrenfest equation (4) 
which corresponds to the regime without condensate, 



lim n(k,x,t) -> \p 



2 = Vi 2 - 

2 1 1 



On the other hand, lim uj — > k 2 which is different from the Ehrenfest cxpres- 

p^O 

sion uj = k 2 + U. Thus, one cannot recover the non-condensate (Ehrenfest) 
description by just taking the limit of zero condensate amplitude in (10), (11) 
and (14). However, one can easily write a unified WKB description which will 
be valid with or without condensate by simply adding U + pto the frequency 
(14). Indeed, for strong condensate U + i o=const and, therefore, it does not 
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alter the ray equations (which contain only derivatives of ui). On the other 
hand, such an addition allows us to obtain the correct expression 

uj = k 2 + U, 

in the limit p — > 0. Summarizing, we write the following equations of the linear 
WKB theory which are valid with or without the presence of a condensate, 

D t n(x, k, t) = 0, (17) 

where 

(18) 

is the waveaction and 

D t = d t + x-V + k-d k , (19) 

is the full time derivative along trajectories and 

x = dkOJ, k = -Vw, (20) 

are the ray equations with 

u = k^/k 2 + 2g + U + p. (21) 

Formula (21) is an important and nontrivial result which can be obtained 
neither from existing general facts about the WBK formalism nor from the 
linear theory of homogeneous systems. 

4 Weakly nonlinear GP equation 

The derivation for the description of the nonuniform turbulence found in a 
BEC system consists of a amalgamation of a WKB method, for the descrip- 
tion of the linear dynamics, and a standard weak turbulence theory (see e.g. 



. 1 ujp 
n(k,x,t) = --p- 



Q q 

UJ 
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[13]), with the noted modification that Gabor transforms are used instead of 
Fourier ones. We will now demonstrate the general ideas of such a derivation 
for the simple case of system where no condensate is present. 
Consider the Gabor transformation of (1): 

idt4 + A$-\^-Uj> + i{V x U)V k $ = 0. (22) 

To calculate the IV'PV' term let us first separate the Gabor transform into its 
correspondingly fast and slow spatial parts, 

i>(x,k,t) = a(x,k,t)e^. (23) 

slow f &s t 
Now by using the inverse Gabor transform 

g(x,t) = J g(x,k,t)dk, (24) 

we find 

I^F^ = e lk x j f(x - x ) e i*o-(fe3+fe2-fei-fe) 

x a*(fci, x )a(k 2} x )a(k 3 , x ) dx dkidk 2 dk 3 . 

(25) 

Note that the slow amplitudes a do not change much over the characteristic 
width of the function / and hence their argument x can be replaced by x. 
Therefore, we can approximate (25) by 



M-< - - — / / • /.• '. k 2 k k • 



e" 

(27) 3d / 2 

x a*(fci, x)a(k 2 , x)a(k 3} x) dkidk 2 dk s . 



(26) 



Here F(k) is the Fourier transform of f(x). Note that for the spatially ho- 
mogeneous systems, e* — > 0, F(k) is just a delta function, 

Inn F(fc) -» S(k). 
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After dropping terms proportional to Ao, equation (22) then becomes 

d t a(k, x) = —2k ■ Va(fc, x) 

-i(k 2 + k- (V x ))a(fe,x)- (V x U)(V k a(k,x)) 

— J F(fc 3 + fc 2 — ki — k) a*(ki 7 x)a(k 2 , x)a(k 3} x) dkidk 2 dk 3 . 

(27) 

This is the master equation formulating the nonlinear dynamics in terms of 
the Gabor amplitudes. This can serve as a starting point for the statistical 
averaging which in turn leads to the weak turbulence formalism. Note that 
this equation can be written in Hamiltonian form, 

with a Hamiltonian function 

H = J {uJk,x - x-V x w/s jX )|afe iX | 2 

+ 2^ 'k^k,x)(ak,x^ 'xa%,x - a l, x ^xak,x) dkdx 
+ J F(k 3 + k 2 -k 1 -k) 

a*(ki, x)a(k 2 , x)a(k 3 , x)a(k 4 , x) dkxdk 2 dk' i dk il 

(29) 

where cj kyX — k 2 + U(x). In fact, such a Hamiltonian description can be 
derived directly, in terms of the Gabor amplitudes, from the Hamiltonian 
formulation of the original GP equation (see Appendix B). 

If a condensate is present in the system, one can also re- write the equations 
in a Hamiltonian form with an identical quadratic part. That is, with a being 
replaced by the normal amplitude, and iv by the frequency of waves, found 
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in the presence of the condensate. It appears that the quadratic part of the 
Hamiltonian (29) is generic in the WKB context. Indeed, let us consider a 
typical Hamiltonian for linear waves in weakly inhomogcncous media [32] 
expressed in terms of Fourier amplitudes a qi and a* 1 

H = J Q (qi , q) a qi a* 1 dqdqi , (30) 

with a hermitian kernel I2(q l7 q) = i?(q, qi) which is strongly peaked at 
q — c\x =0. As we will show in a separate paper [26], this Hamiltonian can 
be represented in terms of the Gabor transforms as 



H = j \oJk,x - X-V x Wfe,a!)|afe,x| 



1 

1 

+ kUk, x )(a kyX V x a* k x - a* k x V x a k:X ) dkdx, (31) 

where a kx are the Gabor coefficients, and LOkx is the position dependent 
frequency, related to J7(q, qi) via 



Wfe, i 



J e- 2lq x n(k,k + 2q)dq. (32) 



Actually, such an expression is a canonical form, even for a much broader class 
of Hamiltonians that correspond to a significant class of linear equations with 
coordinate dependent coefficients [26]. That is, 

H = y*[yl(qi,q)a qi a q S(qi,q)a qi a_ q + c.c] dqdq l7 (33) 

where functions A and B peaked at q — qi =0. 

5 Weak turbulence for inhomogeneous systems 

Now, by analogy with homogeneous weak turbulence, we define the waveac- 
tion spectrum as 

n Kx - (\a(k,x)\ 2 )/F(0), 
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where averaging is performed over the random initial phases. Note that this 
definition is slightly different to the usual definition of the turbulence spec- 
trum in homogeneous turbulence, i.e. the definition constructed from Fourier 
transforms, rikS(k — k') = (a(fe)a(fe')). Indeed, a Gabor transform can be 
viewed as a finite-box Fourier transform, where k = k' in the definition of 
the spectrum and one replaces S(k — k') with the box volume F(Q). 

Multiplying (27) by a*(k,x) and combining the resulting equation with 
its complex conjugate, we get a generalization of (4): 

D t n k ,x = -23 J F(k + fei - fe 2 - fc 3 ) 

x (a*(k, x)a*(k\,x)a(k2, x)a(k 3} x)) dk 1 dk2dk 3 , 

(34) 

with D t = d t +x- V + A; • <9fc . Note, that in the case of homogeneous turbulence, 
using the random phase assumption, in the above equation, would lead to 
the RHS becoming zero. This means that the nontrivial kinetic equation 
appears only in higher orders of the nonlinearity. For the inhomogeneous 
case, the nontrivial effect of the nonlinearity appears even at this (second) 
order. This can be seen via a frequency correction which, in turn, modifies 
the wave trajectories. This effect was considered by Zakharov ct al [31] and 
it is especially important in systems where such frequency corrections result 
in modulational instabilities followed by collapsing events. In our case the 
nonlinearity is "defocusing" and, therefore, such an effect is less important. 
Indeed, in what follows we will neglect this effect as, at sufficiently small 
ratios of the inhomogeneity and turbulence intensity parameters, e < ^, 
wave collision events are a far more dominant process. 
Let us introduce notations 

J k2k 3 = (a*(k,x)a*(k 1 ,x)a(k 2 ,x)a(k 3 ,x)}, 
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and 

1^1% = (a*(k,x)a*(k 1 ,x)a*(k 2 ,x)a(k 3 ,x)a(k 4 ,x)a(k 5 ,x)}. 
Then, we have the following equation for the 4th-order moment, 

A4% fe f = i(£> fci + u k , 2 - w fc , - u K )I k k £? 

+ /(CSP( fe i + fe 1 - fe 2 - fe 3) 

+< i 55 F ( fc 2+fci-fc2-fe 3 ) 

-4lS F ( fe 3 + fel-fc2-fe 3 ) 

-4t£t li? ( fe 4 + fei - A; 2 - fc 3 ) )dk 1 dk 2 dk 3 , 

(35) 

where we denote cDfc = fc 2 + (fc • V^f/). Note that the first two terms on the 
RHS of this equation can be obtained one from another by exchanging k[ 
and fc' 2 , whereas the last two terms - by exchanging k' 3 and k' 4 . To solve this 
equation, one can use the random phase assumption which is standard for 
the derivation of a weak homogeneous turbulence theory and which allows 
one to express the 6th-order moment in terms of the 2nd-order correlators. 
For homogeneous turbulence, the validity of this assumption was examined 
by Newell et al [18,33] who showed that initially Gaussian turbulence (char- 
acterized by random independent phases) remains Gaussian for the energy 
cascade range whereas in the particle cascade range deviations from Gaus- 
sianity grow toward low k values. However, these deviations remain small over 
a large range of k for small initial amplitudes and the random phase assump- 
tion can be used for these scales. Note that the deviations from Gaussianity 
at low k correspond to the physical process of building a coherent condensate 
state. The results of [18,33] obtained for homogeneous GP turbulence will 
hold for trapped turbulence too because inhomogeneity has a neutral effect 
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on the phase correlations. Indeed, according to the linear WKB equations 
the phases propagate unchanged along the rays. Thus we write 



r 123 
J 456 



n in2 n 3 (Fi(FiFZ + F*F*)+ Ff (i^ 1 + F%Fq) 

+ F!{FlFl + FlFl)), 



(36) 



here we have used the shorthand notations, F 2 J = F(0) 5(k\ — fc 2 ) and J||| = 
^klktkl ■ Using this expression in (35) we have 



D 

Dt' 

+ 2 (n fe3 n ki (n kl + n k2 ) - n kl n k . 2 (n k3 + n k J) . 



Notice that the u k terms get replaced by co k , since the (fc • V X U) terms drop 
out on the resonant manifold. Let us integrate this equation over the period 
T which is less than both the slow WKB time 1/e and the nonlinear time 
l/cr 4 . Then, one can ignore the time dependence in n k on the RHS of the 
above equation and we can take k = —X7U = const on the LHS. 

The resulting equation can be easily integrated along the characteristics 
(rays) which in the limit loT — > oo gives 

T klkl = -2{n k3 n ki (n kl n k2 ) - n kl n k2 (n k3 + n ki )] 

5{uj kl + u) k2 - uj k;i - WfeJ. (37) 

Note that to derive a similar expression in the theory of homogeneous weak 
turbulence one usually introduces an artificial "dissipation" to circumvent 
the pole and to get the correct sign in front of the delta function (see e.g. 
[14]). The roots of this problem can be found even at the level of the linear 
dynamics, where the use of Laplace (rather than Fourier) transforms pro- 
vides a mathematical justification for the introduction of such a dissipation. 
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However, in our case there is no need for us to introduce such a dissipation 
because inhomogeneity removes the degeneracy in the system. Substituting 
(37) into (34) we get the main equation describing weak turbulence, the four- 
wave kinetic equation 

D t n k = -[ n k nin 2 n 3 ( — + — — ) S (k + fei - k 2 - k 3 ) 

tv J \n k m n 2 n 3 J 

S (wfc(x) + wi(x) - lu 2 (x) - W3(x)) dkidk 2 dk 3l 

(38) 

where, 

D t = d t + x ■ V + k ■ dk, x = 9few, k = —Voj. 

We can see that the main difference between the kinetic equation for inho- 
mogeneous media and homogeneous turbulence [11-13,22] is that the partial 
time derivative on the LHS is replaced by the full time derivative along the 
rays. Further, the frequency u> and spectrum n are now functions not only of 
the wavenumber but also of the coordinate. 

The same is true for the case when the ground state condensate is im- 
portant for the wave dynamics [13]. The main interaction mechanism now 
become three wave interactions, with the kinetic equation 

D t n = 7r J iVfefcifcal 2 /fei2 <5k-k!-k 2 <L k - Wkl -w k2 dki<2k 2 

-2ttJ \V klkk2 \ 2 h dkidk 2 , (39) 

where f k \ 2 = n^n^ —n^n^ +rik 2 ) • Here, n k , D t and u) are given by expres- 
sions (18), (19) and (21) respectively and the expression for the interaction 
coefficient V kklk2 can be found in [13]. Three- wave interactions always domi- 
nate over the four- wave process when p <~ k 2 (because k <~ 1 and n <C 1). In 
the case /)<P, the relative importance of the three- wave and the four- wave 
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processes can be established by comparing the characteristic times associated 
with these processes. The characteristic time of the three wave interactions 
for p <C k 2 is 

T3w = k 2 ~ d /gn. 

Thus, the 3- wave process will dominate the 4- wave one if the condensate is 
stronger than the waves, i.e. if g > nk d ~ 2 . 

6 Summary 

In this paper, we developed a theory of weak inhomogeneous wave turbulence 
for BEC systems. We started with the GP equation and derived a statistical 
theory for the BEC kinetics which, in particular, describes states which are 
very far from the thermodynamic equilibrium. Such nonequilibrium states 
take the form of wave turbulence which is essentially inhomogeneous due to 
the fact that the BEC is trapped by an external field. There are two main 
new results in this paper. First of all, we have described the effect of the inho- 
mogeneous ground state on the linear wave dynamics and, in particular, we 
have shown that such an effect cannot be modeled by renormalizing the trap- 
ping potential as it was previously suggested in literature. This was done by 
deriving a consistent WKB theory based on the scale separation between the 
ground state and the waves. Our results show that the condensate "mildly" 
pushes the wave turbulence away from the center but it can never reflect it 
(as an external potential would). Note that we established this result only 
for the limit of large occupation numbers described by the GP equation and 
this, in principle, does not rule out a possibility that the the renormalizcd 
potential approach can still be valid in the opposite limit of small occupation 
numbers. Secondly, we showed that the kinetic equation for trapped waves 
generalizes, and one can combine the linear WKB theory and the theory of 
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homogeneous weak turbulence in a straightforward manner. Namely, the par- 
tial time derivative on the LHS of the kinetic equation is replaced by the full 
time derivative along the wave rays, while the frequency and the spectrum 
on the RHS now become functions of coordinate. A suitable definition for 
the coordinate dependent spectrum is given by using the Gabor transforms 
instead of Fourier transforms. It is important to notice that the coordinate 
dependence of the wave frequency has a profound effect on the nonlinear 
dynamics. The resonant wave interactions can now take place only over a 
limited range of wave trajectories which makes such interactions similar to 
the collision of discrete particles. 



Similarly to the case of homogeneous turbulence considered in [13], the 
presence of a condensate changes the resonant wave 



interactions from four-wave to three-wave if the condensate intensity ex- 
ceeds that of the waves. A distinct feature of the inhomogencous turbulence 
trapped by a potential is that if the three-wave regime is dominant in the 
center of the potential well, it is likely to be suddenly replaced by a four-wave 
dynamics when one moves out of the center beyond the condensate reflection 
points where the condensate intensity is decaying exponentially fast. Thus 
the same wavepacket can alternate between three-wave and four-wave in- 
teractions, with other wavepackets, as it travels back and forth between its 
reflection points in the potential well. (The wavepacket reflection points being 
further away from the center than the condensate's own reflection points). 
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Appendix A: derivation WKB equations in presence of 
a condensate 

Let us split <f> into its real and imaginary parts a — $t<j> and b — 3</>. Then 
the equation (9) splits into two coupled equations 

d t a + Ab + 2v ■ Va + -j- ■ V6 + p(2ab + b(a 2 + b 2 )) = 0, (40) 

d t b - Aa + 2ag + 2v ■ V6 - — • Va + g(3a 2 + b 2 + a(a 2 + b 2 )) = 0, (41) 

Q 

where we have used the fact that ^j 2 - = ^ + which follows from (6). 

Gabor transforming our two coupled equations (40) and (41) and using 
Taylor series to represent large-scale quantities, 

g(x ) = q{x) + {x - x) ■ Vg{x) + 0(e 2 ), 

we find 

d t a + Ab + — ^ • V6 + v ■ Va + Q [p(2ab + b(a 2 + b 2 ))] = 0, (42) 
Via 

d t b - Aa • Va + 2v ■ Vb + 2ga + 2iVg ■ d k a 

g 

+ G [{g{'ia 2 +b 2 + a(a 2 + b 2 )))] = 0. (43) 

Where G[f(x)] is the Gabor transform of f(x). We have kept only 0(e) terms 
and neglected the 0(e 2 ) and higher order terms. For generality, we have kept 
the nonlinear term. 



The e° order - As in all WKB based theories we first derive a linear 
dispersion relationship from the lowest order terms. At zeroth order in e, the 
spatial derivative of a Gabor transform is Va = ika which is similar to the 
corresponding rule in Fourier calculus. Then, at the lowest order, equations 
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(42) and (43) become 

d t a - k 2 b = 0, (44) 

d t b + k 2 a + 2ga = 0. (45) 

These two linear coupled equations make up an eigenvalue problem. Diago- 
nalizing these equations we obtain 

dtX = +iuj\, dtH = —iui\x. 

Correspondingly, we find the eigenvectors 

or, re-arranging for a and b 

a = \ + fj,, S=p(A-/i). (47) 

The eigenvalues are given by the dispersion relationship, 

lu 2 = k 2 (k 2 + 2g), (48) 

which is identical to the famous Bogoliubov form [21] which was also obtained 
for waves on a homogeneous condensate in the weak turbulence context in 
[13]. 

Therefore, at the zeroth order, we see that A rotates with frequency — ui 
and fi rotates at +ui. Note that the A and fi are related via 

X*(k)=fi(-k). (49) 

The e 1 order - Let us split the wave amplitudes into fastly and slowly 
varying parts, 

X(x,k,t) = A{x 1 k 1 t)c lk - X+W \ ii(x,k,t) = M(x,k,ty k - x - luJ \ (50) 
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or, in shorthand notation, 



X = Ae+, n = Me-, (51) 

where 

The e + and e~ represent the fastly oscillating parts of the Gabor transforms. 
From (49) it follows that 

A*(x,-k,t)=M(x,k,t). (52) 

Obviously, 

d t \ = + c + d t A, dtn = — iu^ + c~<9 t M, (53) 



VA = ik\ + c + WA + itXVoj, V/z = ikfi + e"VM - iizzVw, (54) 

AA = -k 2 X + 2ic + k ■ VA + e+A/1 - 2t\k ■ Vlu, (55) 
A/z = -fc 2 /z + 2ic + k ■ VM + c+AM + 2tfik ■ Vw, 

<9/cA = c + dkA + ixc + A + itc + Adk^ 1 (56) 
dkft = e~dkM + ixe~M — ite~MdkUJ- 

Our aim now is to derive equations for d t X and d t \i- However, due to the 
relationship (49) it is sufficient to derive an equation for only one of the two, 
for example A. From (46) we find 
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After substituting our equations for d t a and d t b, (42) and (43), and making 
use of the relationships (47) the equation for A acquires the following form: 



<9 t A = A 
+ AA 

+ A/i 



iV g ■ Vuj ik 2 g 
2k 2 g + ~o7~ 
iuj 



2uj 



+ VA- 



-Vg ■ dkX 



iVuj iujVg ik 2 Vg 
— 2v — 



k 2 



2k 2 g 



2ujg 



iVg ■ Vuj ik 2 g 
2k 2 g + ~o7~ 



Vp- 



iVuj iujVg ik 2 Vg 



k 2 2k 2 g 2ujg 



IUJ 

2k 2 



2uj 



-Vg ■ d k p - ML. 



Here the nonlinear term ML is given by 

ik 2 

ML = Q [p(2ab + b(a 2 + b 2 ))] - [(g(3a 2 + b 2 + a(a 2 + b 2 )))] . 

2oJk 

Note that we have neglected Co in the above expressions because, according to 
the dispersion relationship (48), it is of the order of p which is 0(e 2 ) by virtue 
of (7). We will also drop the nonlinear term in the subsequent calculation. 

Our next step is to eliminate the fast oscillations associated with the 
Gabor transforms and derive an equation for \A\ 2 . This in turn will lead to 
a natural waveaction quantity which can be used to describe the behavior of 
our wavepackets in phase space. Using (53-56) we obtain 



^ , A «Vg • Vw ik 2 g 
dtA = A [—2k^g- + — 

+ [ikA + VA + itAVuj] ■ 



IUJ 

iVuj iujVg 



k 2 



2v 



2k 2 g 

iuj ik 2 

~~2~k? ~ 2oj~ 

k 2 ^ „ A ik 2 A „ itk 2 A^ „ 
Vg ■ OkA x ■ Vg Vg ■ OkUJ- 

UJ UJ UJ 



ik 2 Vg 
2ujg 



+ [-k 2 A + 2ik ■ VA - 2tAk ■ Vuj] 



(57) 



Please note that all the terms involving M drop out. This stems from the 
fact that, in deriving an equation for A, we have had to divide through by 
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o + . Therefore, any terms involving M will result in a factor 



e /e + = e 



-2iut 



Thus, after time averaging over a few wave periods, all the M terms drop 
out. 

Expanding out equation (57) we find the 0(1) terms cancel out and using 
the dispersion relationship (48) wc find 

d t A = d k uj ■ VA + ^k -Vq-Voj- d k A + U, (58) 

where 

tAw, „ tk 2 A, „ „„, k 2 A „ tk 2 A^ n 
J = -r^-fe • Vcj H k ■ voj — 2Ak ■ v x ■ vg \7g ■ Ok^ 7 

At this point let us drop the nonlinear term and concentrate on the linear dy- 
namics. Multiplying (58) by A* and combining it with the complex conjugate 
equation the J terms cancel, leading to 

d t \A\ 2 - d k cu ■ X7\A\ 2 + Vuj ■ d k \A\ 2 = ^j^k ■ Vg. (59) 

A similar equation for \M\ 2 can be easily obtained by replacing k — > — k in 
(59) and using (52), 

8 t \M\ 2 + d k tu ■ V|M| 2 - • 8 k \M\ 2 = -^f^ ' Vg. (60) 

The LHS of this equation is the full time derivative of \M\ 2 along trajectories. 
If \M\ 2 were to be a correct phase-space waveaction, the right hand side of 
this equation would be zero, however, this is not the case. We find the correct 
waveaction n(x, k, t) by setting 

\M\ 2 = a(x,k)n(x,k,t), 
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and finding such a(x, k) that the the full time derivative of n(x, k, t) is zero. 
This leads to the following condition on a, 

dkUJ ■ Va — Vw • dkoi H p-fe • V g = 0. 

By choosing a = k x g Y and substituting it to (6) we find x = 2, y = — 1. 
Therefore the correct form of the waveaction is n = p-|M| 2 . Summarizing, we 
have got the following transport equation for the waveaction n in the linear 
approximation, 

D t n(x, k, t) = 0, (61) 

where 

D t = d t +x-V + k-d k , (62) 
is the full time derivative along trajectories and 

x = dkOJ, k = -Vw, (63) 

are the ray equations with 

lo = k^/k 2 + 2g. (64) 

Obviously, the dynamics in this case cannot be reduced to the Ehrenfest 
theorem with any shape of potential U. Therefore, approaches that model 
the condensate effect by introducing a renormalizcd potential are misleading. 

Finally, it is useful to express the waveaction n in terms of the original 
variables, 

2 

(65) 

It is interesting that such a waveaction is in agreement with that found in [13]. 
In fact in [13] the homogeneous case with non-zero nonlinearity (e = 0, a ^ 0) 
was considered. This is the opposite limit to the one we have considered above 
(where e ^ 0, u = 0). 



n(k, x, t) 



1 ujp 



ik 2 



30 Lvov, Nazarenko, West 

Appendix B: Hamiltonian formalism for spatially 
inhomogeneous weak turbulence. 

Let us start with the GP equation written in the Hamiltonian form: 

l df* = (66) 

The Hamiltonian for the GP equation (1) coincides with the total energy of 
the system: 



n= I '/r( W x | 2 + i|!P x | 4 + ^(x)|!? x | 2 ). (07) 



Let us first consider the case without a condensate. Applying the Gabor 
transformation to (66) we get 



But if wc notice that 



d f 8H 



6V X 6& 



we obtain 



x,k 

at SKm 

Thus, the time evolution of the Gabor transformed quantity is governed by 
the Gabor transformed Hamiltonian equation. However, we would like to 
obtain the equation of motion in Hamiltonian form without the Gabor trans- 
formation. Let us re-write (69) in terms of the slow amplitudes a defined in 
(23) 

^a k , x = //(x-x0^x'. (70) 
Now, let us express the Hamiltonian (67) in terms of the slow variables a, 

J eW 1 -^* (-a kl , x (-fc 2 2 + 2*fc 2 V) a £ 2!X + (7(x) akl ^ 2!X ) drdk^ 

+ J e i ( fcl+fa - fe -^ )x Ok 1 ,xOk a ,xa^ >x oJ 4iX drdkidkadksdk 4 . 

(71) 
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Here we have integrated by parts IV'Z'I 2 and, while calculating the Laplacian 

of & in terms of slow variables, have kept only the first order gradients in 

a k , x - Substituting (71) into (70) allows us to re-write this equation as 

d 5H 
l 8i a ^ = to£> (?2) 

where the filtered Hamiltonian H can be represented as 

H = J f{x-x')e i{kl - k2)x ' 

x (— ciki,x'(^2 + 2«fc2V)a£ 2 . x + y(x')aki,x'aik 2i x) drdr'dk 1 dk 2 
+ J F(k 1 + k 2 k 3 k 4 ) a kl>x fl k2jX a^X 4 , x dxdk 1 dk 2 dk 3 dk 4 , 

(73) 

and F(k) is the Fourier transform of the f(x). 

Expanding U(x') as U(x) + (x — x')VU(x) and taking into account that 
(x — x') can be interpreted as —id k2 e lk2 ^ x ~ x \ we have 

H = J ((fc 2 + C/( a; ))| afej;c | 2 _l(VC/(x))# fei;c 9 fc ^ ;E +ifc afe!;c V4 !;c 
+ l -{VU{x))$l x d k 4> k , x - <h„. X «•,..,) dkd-K 
+ J F(ki + k 2 k 3 - k 4 ) ak 1; xak 2 ,xak 3; x a k4,x d^dk ± dk 2 dk z dk Al 

Since oo k ,x = k 2 + U (x) we can represent the above formula as 

H = j {{uJk,x - xVuj k , x )\a k , x \ 2 + ^(V x u; fex ) (a* ktX Vka.k,x ~ a ktX V k a* kx ) 
+ ^(^k^k,x) (a k , x V x a* kx - a k , x V k a* k x ) ) dkdx 
+ J F(ki + k 2 - k 3 - k 4 ) a kl xak 2 xak 3 x a k 4 x dxdk 1 dk 2 dk 3 dk 4 , 

(74) 

Now, we will show that if a condensate is present then the quadratic part 
of the Hamiltonian can also be written in the same canonical form as in (74) . 
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Let us start from the equation (58) for A 



d t A = d k uj ■ VA + -^fc ■ S/g - Vw • d k A + iJ, (75) 



with 



tAu, „ tk 2 A, „ „„, k 2 A „ ifc 2 ^l„ „ 

J = —rr^-k ■ \7u> H fe • Vw — 2/lfc • t; a; • Vg vg ■ OkUi, 

k z ui ui ui 

Expression (65) for the waveaction in this case allows us to guess the form of 
the normal variable, 



_ ^g^k.x iLUh w 

Note that this expression is consistent with the waveaction considered above 
for the case with no condensate. This can be checked by taking the limit 
p — > 0. In terms of normal variable afe jCC equations (58) and (75) acquire the 
following form: 

<jfe,x = iuk, x ak,x + 9fcWfc,xVa fc;X + Vujk, x dkak, x 
- 2iau, x k ■ v - ia ktX x ■ Vw fc ^. 

(76) 

This equation can be represented in the form of a Hamiltonian equation 
of motion with a quadratic Hamiltonian as in (74) when the frequency is 
replaced by its Doppler shifted value, 

ui — > ui + 2k ■ v. 

Note that the Doppler shift does not enter into the equation for the wave- 
action because it leads to terms that are of second order in e and therefore 
should be neglected. 
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